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Summary. We propose an iterative finite element method for solving non-linear 
hydromagnetic and steady Euler's equations. Some three-dimensional computational 
tests are given to confirm the convergence and the high efficiency of the method. 



1 Introduction. Statement of the problem 

The understanding of plasma equilibria is one of the most important problems 
in magnetohydrodynamics and arises in several fields including solar physics 
and thermonuclear fusion. Such an equilibria is often governed by the well 
known steady hydromagnetic equations 



which describe the balance of the Lorentz force by pressure. Here B and p are 
respectively the magnetic field and the pressure. 

Notice that system (l)+(2) is quite similar to steady inviscid fluid equa- 
tions 



curlB xB+Vp = 0, 
div B = 0, 



(1) 

(2) 



v.Vv + Vp = 0, 
div v = 0. 



(3) 
(4) 



This analogy is due to the vectorial identity 
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v.Vv — V— — = curlv x v. 

System of equations (l)+(2) must be completed with some boundary condi- 
tions on B and p. Physical considerations suggest to prescribe the boundary 
normal field component: 

B.n = g on dfi (5) 
where g satisfies the compatibility condition / g = due to the equation 



divB = 0. Defining the inflow boundary as r~ = {x € fi, B(x).n(x) < 0}, 
one can also prescribe the normal component curl B.n of the current density 
and the pressure pon 

curl B.n = h on r~ , (6) 

p = po on r~. (7) 
One can notice that if the pressure is neglected, equations (l)+(2) become 

curlBxB = 0, (8) 
divB = 0. (9) 

Equation (8) means that the magnetic field and its curl, which represents the 
current density, are everywhere aligned. The magnetic field is said Beltrami 
or force-free (FF). A usual way to tackle the problem (8) + (9) consists to 
rewrite equation (8) into the form 

curlB = A(x)B, (10) 

where A(x) is a scalar function which can be a constant function or can depend 
on x. In the former, the B field is said linear FF. In the latter, it is said non 
linear. 

Some partial results concerning existence of 3D solutions of equations (l)+(2) 
in bounded domains are given in [1] and [11]. Linear force- free-fields were 
studied in [4] . For the existence of non-linear ones the reader can refer to [5] , 
[3]. 

The numerical solving of equations (l)+(2) and equations (8) + (9) is of 
importance in magnetohydrodynamics studies and in solar physics. As it is 
known, the reconstruction of the coronal magnetic field has is of a great utility 
in observational and theoretical studies of the magnetic structures in the solar 
atmosphere. In this paper, we propose an iterative process for solving these 
equations (section 2). A finite element method is proposed for solving each 
one of the arising problems. 



2 An iterative method for the magnetostatic system 



Our objective here is to expose an iterative method for solving the non-linear 
equations (l)+(2) in a bounded and simply-connected domain. The starting 
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idea of the method consists to split the current density uj = curl B into the 
sum 

U) = W|| + W_L, (11) 

where the vector field u;|| = m( x )B is collinear to B, while u>j_ is perpendicular 
to B. The problem is decomposed formally into a curl-div system on B(x) 
and two first order hyperbolic equations on /u(x) and p(x). 
More precisely, writing W||(x) = ^i(x)B(x) where [i is a scalar function and 
taking the divergence of (11), gives 

B.V/i = -divu;^. (12) 

Notice that the pressure satisfies a similar equation since 

B.Vp = 0. (13) 

Equation (1) becomes 

x B = -Vf>, (14) 

which means that w_l(x) = — }■ Vp(x) x B(x) if |B(x)| ^ 0. 

In consideration of these remarks, we are going now to propose an iterative 
process to solve non-linear systems (l)+(2). In this process the transport 
equation (12) is perturbed by adding an artificial reaction term. Namely, we 
construct a sequence (B( n ',p( n ') n > as follows: 

• The starting guess Bo € if 1 (J?) is chosen as the irrotational field associ- 
ated to g defined by 

curlB = in J?, divB = in fi and B .n = g on dfl. (15) 

This is a usual problem which can be reduced to a scalar Neumann problem 
since the domain is simply-connected. 

• For all n > 0, p^ is solution of the system 

BKVpW + VP (n) = VP (n ~ 1) in V, 

= po on on, ( ' 

where r/ is a small parameter and p^ -1 ' = 0. 

For all n > 0, u> ± ^ = , ,\ ln Vp (n) x and w M ( n ) = m (,i) B(™\ where 
- - 1 - |B(")| 2 11 

satisfies 

' B(").VAi (n) + eAi (n) - -divw_L (n) + eM (n_1) in ^, / 17l 
Here ^i^ 1 ) = and e is a small parameter. 
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• For all n > 0, B(™ +1 ) = B + b( ,l+1 \ with b<" +1 ) solution of 

curlb(" +1 ) = + Vq (n) in Q, 
divb(™ +1 ) = in Q, 

b (n+i) _ n = o on an, 

where u/ n ) = + u±^ while is solution of the Laplace problem 

-Aq {n) = divu} (n) in [2, and q {n) = on dfl. (18) 

Notice that the appearance of the correction term Vq^ is due to the fact 
that div (a/™)) is not zero in general. 

The convergence of this iterative process is not an easy matter. We conjecture 
that it converges if h is sufficiently small and |B (x)| > c > in Q for some 
constant c > 0. Nevertheless, in the case of linear force-free fields (in that case 
the algorithm is simplified since at each iteration p^ — 0, u>±^ — and 
/l/™) is a fixed real) Boulmezaoud and Amari [6] proved that this process is 
super-convergent. The proof of convergence in the general case is not given 
and remains an open question. 

Notice that the same algorithm can be used for computing linear or non- linear 
force- free fields which are solutions of (9)+(10), provided that the computation 
of the pressure and the vector field uj±^ are dropped. 



3 Finite element discretization 



Here we give a short description of the finite elements methods we use for 
solving problems arising in the iterative process exposed above. Observe first 
that at each iteration of the algorithm one should solve two problems: 
(a) A reaction-convection problem of the form: find u solution of 



div (uB) + au = f in Q, 

u = h on r~ 



(19) 



(b) A vector potential problem: find the pair (b, q) satisfying 



curlb = j in Q 1 

div b = in i?, 

b.n = on df2, 

q = on dfl. 



(20) 



We begin with the approximation of (19). 

It is well known that the direct application of a Galerkin finite elements 
method to the singularly perturbed problem (19) may lead to the appearance 
of spurious oscillations and instabilities. In the two last decades, several meth- 
ods were proposed to remove this drawback (especially in the two dimensional 
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case). Among these methods, one can recall the streamline diffusion method 
(sec Brookes and Hughes [8], see also, e. g., Johnson et al. [10]), the discon- 
tinuous Galerkin method (see Lesaint [12]) and bubble functions methods (see, 
e. g., Brezzi et al. [7]). Here we shall use the method of streamline diffusion. 

Thus, let us consider a family of regular triangulations (^h) of fi. The 
discrete problem we consider is 



Find Uh <E Wh such that 
a h (u h ,w h ) =£h(w h ), Vw h G W h , 



where 

a h (u h ,w h ) = / (B.Vu h + o-Uh).(wh+5hB.Vw h )dx- / u h w h (B.n)dx, 

JO Jr- 

h(wh) = / /(xXB.Vwfc + S h w h ) - / a a w h (B.n)dx. 

Jn Jr- 

Here Wh stands for the finite elements space 

W h = {v h € H\Q); v\ K e P k (K), VK E ST h }, 

where for each K e 3?h, ^k{K) denotes the space of polynomials of degree 
less or equal k. 

One can prove that the problem {S?h) has a unique solution u h € Wh when 
5h<y < 1. Moreover, if Sh = ch for some constant c and if B e L°°(J?) 3 n 
H{div ; Q) and u e H e+1 ([2) for some I > 1, then 

(l-S h a)\\\u-u h \\\ <Ch l+1 ' 2 \\u\\ H ^ (n) , (21) 

where |||w|||^ = 5h\\B.Vw\\ 2 L2{{2) + a\\w\\ 2 L2{n) + || IB.nl 1 / 2 ^!!^^. 
Now, we deal with the approximation of the curl-div system (20), which can 
be dispatched into two problems: a variational problem (=2) in terms of b and 
the fictitious unknown 9 = 0, and Laplace equation (18) in terms of q. We only 
deal with the approximation of b, since we shall see that the computation of 
the q is useless. Denote by iJ(curl ; Q) the space 

#(curl; Q) = {v e L 2 (f2) 3 ; curlv e L 2 (Qf} 

equipped with its usual norm. The statement of problem (20) suggests the use 
of an iJ(curl ; Q) approximation. Define M the space 



M = {v e H^Q), [ v dx = 0}. 



Let Xh C H (curl ; J?), Mh C M two finite-dimensional subspaces and set 
V h = {v h e X h - (v h , = 0, V/ifc G M h }. 



We make the following assumptions 
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(J4?i) the inclusion {V^h, M/» € Mh} C Xh holds, 
(^2) there exists a constant C such that 

||v h ||o,fi < C||curlv h || , r2 , Vv h e V h . 
The discrete version of problem (20) writes 

' Find (bh,0h) € Xh X smc/i as 
. . ^ Vv h e A", l; y curl b h . cur 1 v h rfa; + ^ v h .\79 h dx = J j.curlv h da;, 



e Af fc , / bh.V// h = 0. 

'12 

According to Amrouche and al. [2], the problem (=S/() has one and only one 
solution (hh,6h) with 8 h — 0, and 

!|b-b h || ff(curl:fi) <C inf ||b-v h || ff(curl . r2) . (22) 

A simple manner for constructing the spaces Xh and Mh is to use the H (curl ) 
conforming elements of Nedelec [13] (see Amrouche and al.). In that case, the 
following estimate holds 

||b-b h || ff(cnr i ;fl) <CT/{|b|^ 2 + |b|, +1 ,, 2 }, (23) 

which is valid if b e H e+1 (Q). 

An important feature of the discrete system (&h) is that only the discrete 
vector field bh is really unknown. Actually, we know that Oh — . This 
property can be exploited from a practical viewpoint to reduce the discrete 
system to a smaller one by eliminating Oh- In term of matrices, the system 
writes 

■ A curl B T\ / X \ ^ /(jcurV 

B \Y -{ 



(24) 



where A curl is a symmetric and positive square matrix (A curl is not definite 
neither invertible). We can state the following 

Lemma 1. Let A be a square positive, definite and symmetric matrix having 
the same size as A. Then, the pair (X, Y) is solution of (24) if and only if 
Y = and X is solution of 

(A curl + B T AB)X = C curl . (25) 

Remark 1. In Lemma 1, the matrix A curl and the RHS C curl are not arbitrary. 
Indeed, if G denotes the matrix of the operator V : Mh — > Xh, then necessarily 
qT jsfuri _ q an( j qTqcuH _ q These identities are the discrete counterpart 
of the continuous relations div (curl .) = and div j = 0. 

A serious advantage of the new system (25) comparing with (24) is that num- 
ber of unknowns is reduced. 
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4 Computational tests 

In this last section, we expose some computational results we obtain with a 3D 
code. This code use the iterative method and the finite elements discretization 
exposed above to solve problem (f )+(2) and problem (f0)+(9). We compare 
the exact solution and the numerical solution and we show the behavior of 
the errors in terms of h. Two exact solutions are used for the tests. 

• Test 1 (a non-linear force- free- field) . 

Let (r, 9, z) the cylindrical coordinates with respect to a point (a;o,yo>0) 

(x = -3 and y = -3). The pair (B,p) is given B = —i=(e g + e z ), p(x) = 

1 

0. This is a non-linear force- free field with A — — . Table 1 shows the be- 

2r 

havior of the residue ||B(™ +1 ) - B'"' \\ ,n and the product curlB<") x B<™) 
versus the iteration number. This example illustrates the superconvergcncc 
of the algorithm. 



Table 1. Evolution of "^"^m and HcurlB< n > x B^"> || 



n 


l| B (-+i)_ B («) | 


k " ||curlB (n) x B^Hoc 




l|B<"'|| ,n 







0.09912 


6.740e-15 


1 


0.00566 


0.06781 


2 


0.00036 


0.01939 


3 


2.644e-05 


0.01910 



• Test 2: (Bennet pinch) . The pair (B,p) is given by 

, , A 2A . . ,1 + Xk 2 (x 2 + z 2 ). 
B = \7A x e y and p = -c 2A with A = - ln( ^- -). 

2 ZK 

In table 2, the relative L 2 errors on B and p after convergence of the 
algorithm are shown. These error decreases as ft 1-8 , which confirms the 
high accuracy of the method. 



Table 2. Relative errors on B h and p h in norm L 2 (test 2). 



. I|B-B h lo.o ||p-Ph Wo.n 

n l|B|| ,n ||p||o,fi 

0.69282 0.03837 0.08648 

0.23094 0.00492 0.01102 

0.13856 0.00191 0.00396 

0.09897 0.00108 0.00201 
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Fig. 1. Superposition of the the exact and the numerical solutions in the case of 
test 1 on the left and in a (x — z) plane 2D cut for the test 2 on the right. 
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